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The Hartree-Fock-Bogoliubov approximation is very useful for treating 
both long- and short-range correlations in finite quantum fermion sys- 
tems, but it must be extended in order to describe detailed spectro- 
scopic properties. One problem is the symmetry-breaking character of 
the HFB approximation. We present a general and systematic way to 
restore symmetries and to extend the configuration space using pfaffian 
formulas for the many-body matrix elements. The advantage of those 
formulas is that the sign of the matrix elements is unambiguously deter- 
mined. It is also helpful to extend the space of configurations by con- 
straining the HFB solutions in some way. A powerful method for finding 
these constrained solutions is the gradient method, based on the gener- 
alized Thouless transformation. The gradient method also preserves the 
number parity of the Bogoliubov transformation, which facilitates the 
application of the theory to systems with odd particle number. 



1. Introduction 

Soon after the seminal paper describing the microscopic theory of supercon- 
ductivity by Bardeen- Cooper- Schrieffer (BCS) [T, Bohr et al. [2] found an 
analogy between the excitation spectra of nuclei and those of the supercon- 
ducting metallic state and pointed out the role of pairing correlations in the 
low excitation spectrum of atomic nuclei. As self-bound fermionic systems, 
nuclei are unique in requiring for their theoretical description the inclusion 
of both long- and short-range correlations. The longest range correlations 
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may be treated in the Hartree-Fock (HF) approximation with a suitable 
effective Hamiltonian. In the simplest theory that includes pairing, the 
pairing correlations are introduced through the BCS approximation defin- 
ing the pair amplitudes from the time-reversed orbital wave functions of the 
HF theory. However, in many situations the HF/BCS wave functions are 
not the variational minima in the complete space of wave functions defined 
by the Bogoliubov transformations. For this reason contemporary calcula- 
tions of nuclear structure based on the mean-field approximation (see [5H5] 
for recent reviews) largely follow the Hartrce-Fock-Bogoliubov (HFB) for- 
mulation of theory; see Refs [SI [7] for details in the nuclear physics context. 

The atomic nucleus is a mesoscopic system where the broken symme- 
try implied by the BCS or HFB wave functions is just an artifact of the 
mean field approximation. An improved description of physical properties 
requires techniques beyond the mean field, like particle number symmetry 
restoration or fluctuations in the BCS order parameter. Those techniques 
were developed in the 1960's [8HTT] and applied to a variety of situations in 
nuclear physics [3H1] . Recently, they have been exported to several branches 
of physics [12] and quantum chemistry [13] . Other approaches based on the 
Random Phase approximation and derivatives are also popular (see Y.R. 
Shimizu contribution to this Volume and Ref [14 ) However, technical dif- 
ficulties still remain in its practical implementation, especially in systems 
where time reversal symmetry is broken. One of the difficulties is evaluat- 
ing the sign of matrix elements between two general HFB wave functions. 
The sign is relevant because it determines the interference pattern of those 
linear combinations of mean field wave functions typical of theories for sym- 
metry restoration and/or configuration mixing. The proof that the sign of 
the matrix elements is well defined was given in the past [15j but a general 
and robust methodology to determine it in practice was not available until 
a new method based on pfaffians was introduced [16]. The generalization 
to systems with an odd number of particles (to be denoted odd-A systems) 
has been given recently |17j and our methodology will be discussed below. 

The HFB theory defines a minimization problem that raises the practical 
question of finding the minimum of an energy function that depends on a 
large number of variables. Traditionally the equation for the gradient, ie. 
the derivative of the energy function with respect to all the variables, is 
set equal to zero and the resulting HFB equations are solved iteratively. 
However, it has been long known that there can be severe difficulties with 
this approach, as may be seen in Fig 5.3 of the textbook by Ring and 
Schuck 6 ). The approach using the gradient directly is more stable, and 
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we have taken this path in our group at Madrid to develop efficient codes 
based on a second-order treatment of the gradient. One situation where 
the gradient method has obvious advantages is in treating systems with an 
odd number of particles, discussed in Section|3]below. It is also much easier 
to treat a large number of constraints in the gradient method. This will 
facilitate the extensions of the HFB theory discussed in Section [5] below. 

2. Sign of HFB overlaps with the pfafRan technique 

The problem of calculating the overlap of two HFB wave functions was 
first considered in the 1960's [10] in the context of symmetry restoration. 
The formula derived there involves the square root of the determinant of a 
matrix built with the Bogoliubov amplitudes U and V of the HFB states 
involved. The presence of the square root implies that the sign is undefined. 
However, if time reversal is preserved, Kramers degeneracy implies that the 
determinant in the overlap formula is the square of a number and its sign 
is usually assigned to the overlap (without proof). For general HFB states 
it can be proven |15j that the eigenvalues of the matrix in the argument 
of the determinant are doubly degenerate implying that the determinant is 
again the square of a number. 




Re < w| R 2 (ot)| iv' > 



Fig. 1. Sketch of the real and imaginary parts of a typical overlap of the form 
(w\R z (a)\w') . Filled circles are the values of the overlap; open circles the same but 
with the opposite sign. 



To illustrate the sign problem we present in Fig [2] a sketch of the real 
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and imaginary parts of a typical overlap of the form (w\R z (a)\w') where 
the angle a varies between and 2n. Realistic examples are presented and 
discussed, for instance, in Ref [TB]. In our sketch plot, two sets of points 
are depicted. The filled circles represent the overlaps obtained on a discrete 
mesh of a values. The open circles are the same overlaps but with opposite 
sign. The lines joining the points are plotted to guide the eye. The overlaps 
are used typically in integrals in a (see [19) for examples). From the plot 
it becomes clear that if the procedure to identify the sign is not robust 
(usually arguments based on continuity of the overlap as a function of a 
are used) one can easily jump onto the wrong curve when the modulus of 
the overlap is small. At first sight it could be argued that the error in the 
integral is going to be small as the jump takes place in the region of small 
overlap moduli but continuing in the wrong curve leads to large values of 
the overlaps with the wrong sign. 

An unambiguous evaluation of the sign of the overlap between two HFB 
wave functions was first achieved in Ref |16] . That expression for the overlap 
was derived by the coherent fermion state technique, resulting in a pfaffian 
of a matrix related to the Bogoliubov transformation matrices. While this 
solves the problem for fully paired HFB wave functions, the matrix ex- 
pression can become singular in the HF limit. Other pfaffian expressions 
addressing this and other problems related with the use of different finite 
bases for different states were subsequently found [20] , The limitation in 
these approaches is that only fully paired HFB wave functions are allowed 
and the method is restricted to systems with even number parity. Re- 
cently, a method that uses the expression of the standard Wick theorem 
for mean values of fermion operators in the vacuum in terms of a pfaffian 
has permited the extension of the previous result to odd-A systems [FT] , 
Other treatments of odd-A systems [HI |2T] require the Generalized Wick 
Theorem (GWT) [11] and lead to more elaborated expressions. 

The results obtained in [17] are based on a result for the expectation 
values of fermion operators in the vacuum. The method may be understood 
more easily with an example. If are fermion creation or annihilation 
operators satisfying the standard commutation relations, the standard Wick 
theorem says that the following mean value with respect to the vacuum 

(-\Pl020afi4\-) = ri 2 rzi - r lz r 2i + r 14 r 23 
is given in terms of the contractions = (—\/3i/3j\—). On the other hand, 
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the pfaffian of a general 4x4 (skew-symmetric) matrix is given by 



pf 



/ r i2 ri3 r 14 \ 
-ri2 r 23 ?"24 
-rn -r 2 3 r 34 
\-r u -r 24: -r 3 4 / 



ri2?~34 - ^137-24 + ri 4 r 23 . 



This is exactly the same expression obtained for the above expectation 
value. This suggests the following result: 



(\p 1 ...p P p 1 ...p Q \)=pf(S ij ) 



(1) 



where Sij is the skew symmetric (P + Q) x (P + Q) matrix such that Si 
i < j are all the possible contractions 



(\p k pi\) i,j = l,...,P(k,l = 

(\0 k Pr\) i=l,...,P,j = P- 
(\PrP a \) i.J P ■ 1 /' 



1,...,P) 

-l,...,P + Q(k = l, 
-Q(r,s = l,...,Q) 



,P; 



(2) 
(4) 



We have also introduced another set of fermion operators f3i that are pre- 
sumably related to the /3.; by some canonical transformation. The proof of 
this result can be easily obtained using recursion relations and can also be 
easily extended to finite temperature systems |23j . 

The formula Eq. ([1]) can be readily applied to the problem of computing 
overlaps between two HFB wave functions by noting that such HFB states 
can be written as 



detC 



w 



n 



■Plp2..-hn\) 



(5) 



ct=l 



where the normalization factor in front of the product of quasiparticle an- 
nihilation operators Pi contains the occupancies v a and the determinant of 
the third transformation in the Bloch-Messiah theorem [6l [7] and is con- 
structed to give a normalized \w). An immediate application of this result 
is the formula for the overlap of two HFB states including a canonical trans- 
formation operator 1Z (as the ones that appear when symmetry operations 
are applied to the system) acting on one of the states 



/ /\ / ,„detC*detC" , 



V T U V T R T V* 
-V'^RV WV* 



(6) 



a See |16| for basic results and bibliography concerning pfaffians and 1221 for numerical 
and symbolic techniques. 
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where the matrix R is the representation of the canonical transformation op- 
erator 1Z on the linear Fock space generated by the creation and annihilation 
operator c\ and Cj in some convenient basis, namely IZ.^Hr 1 — J^j Rij c j- 
A general multi-quasiparticle overlap including a canonical transforma- 
tion TZ is easily obtained with the previous formalism |17j 

M^.-.^Kfft ...#t K) = ( _ 1) »(_ 1) .(r-i)/ 2 det(7*det^ ^ (?) 

1 la ^Oi^a 



pf 



V T U V T ^ V T R T q' T V T R T V* 
-p*V q*p f q*R T q' T q*R T V* 
-q'RV -q'Rq< p'q' T p'V* 
-V'^RV -V'tiZqt -y'tp' T t/'ty'* 



(8) 



For this expression to make sense both r and s must have the same num- 
ber parity. The objects p and q (p' and q') are matrices of dimension 
r x 2n (s x 2n) with matrix elements p^ jm = V mtJlj and ? ftm = U mllj . This 
expression has the advantage over the direct application of the general- 
ized Wick's theorm [11] that it avoids the combinatorial explosion of terms 
to be evaluated. Namely, (r + s — 1)!! contractions have to be computed 
if the multi-quasiparticle overlap is evaluated by the generalized Wick's 
theorem. To give an idea of the complexity brought about by the combina- 
torial explosion, let us just mention, for instance, that in the evaluation of 
the Hamiltonian overlap of two quasiparticle excitations built on top of an 
odd-A system, overlaps with ten quasiparticles are required. The number 
of terms to be considered if using the GWT would be 9 !! = 945. If two 
independent two quasiparticle excitations are considered in each isospin 
channel the number of quasiparticle operators increases by four and the 
number of contractions goes up to a whooping 13 !! = 135 135. 

3. Gradient method for the HFB equation of odd-A systems 

Systems with an odd number of particles are less studied from a theoretical 
perspective than even-even systems. Several circumstances could explain 
this imbalance and we now discuss two of them. At the BCS level the wave 
function of an odd-A system is given by [6l [7] 

\<Pko) = a t Q II ( Ul + v i a f a l)\-) 

where the orbital labeled ko is "blocked". As a consequence, this orbital 
acquires an occupancy of one and its time reversed companion fco becomes 
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empty. Another consequence of blocking, the fact that the odd-A BCS 
state is no longer invariant under time reversal, makes it more difficult to 
solve the BCS equations. The Hartree-Fock (HF) and pairing fields also 
acquire time-odd components which must be included in the HFB energies 
and minimization procedures. 

To avoid dealing with the time-reversal breaking issue, people have made 
use of the equal filling approximation (EFA). It amounts to replace the 
density matrix and pairing tensor of a blocked orbital fco by a linear com- 
bination with equal weights of the density matrices and pairing tensors of 
the orbitals fco and fcdB This approximation was widely used even before 
the whole procedure was justified as a variational problem on the energy 
of an statistical admixture of the fco and fco blocked states [H] . Although 
this procedure gives results which are very close to the real blocking when 
the time-odd HF and pairing fields are neglected [25], the differences with 
real blocking can amount to a few hundred KeV and therefore are relevant 
for the determination of spin and parities of the ground and excited states. 

To deal properly with odd-A systems the preferred alternative is the 
HFB approximation with full blocking. The situation becomes even more 
involved than the BCS case because now the odd-A wave function is given 

by 

l</Vo) = alM) 

where a} is the quasiparticle creation operator of the quasiparticle labeled 
/L«o and \4>) is the wave function of an even number parity reference system. 
The reference wave function |</>) is the vacuum of all the quasiparticle an- 
nihilation operators a M , i.e. ctp|</>) = 0. On the other hand, \4>fj, ) is the 
vacuum of the set of quasiparticle operators 

ai,..., a^ _i,a^ ,a Mo+ i, ...,a N . 

The new quasiparticle vacuum can be obtained from the old one [261428] 
by swapping the column hq of U and V. This "swapping" procedure is 
not very easy to incorporate into a practical implementation of the HFB 
method for odd-A systems. This is important from a practical standpoint 
because of odd-A systems typically require many HFB calculations with 
different starting wave functions in order to insure that the ground state is 
reached [29[ 130] , As a consequence, it is very important to have a robust 

b For spherically orbitals the linear combination runs over the 2j + 1 sub-levels with 
weights 1/(2.7 + 1) 
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and efficient method for solving odd-A systems for global applications such 
as the construction of theoretical mass table [25J EU E2] • 

In the context [53] of generalizing the approximate second order gra- 
dient method of [34] it was realized that the "swapping" in the U and V 
amplitudes can be easily incorporated into the formulas. The argument is 
as follows: the most important object in the HFB method is the generalized 
density matrix 



that is given in terms of the unitary Bogoliubov super-matrix 

w=(?,Vl do) 



(11) 



vu\ 

and the generalized quasi-particle density matrix 

(01/3^10) {4>\^f} v \cl>)\ = (0 0' 

When dealing with a blocked HFB state |^ D ) the generalized quasi-particle 
density matrix becomes 

/TTD \ _ / WHO Mo) (^Vo 

where the diagonal matrices Mo and I Mo have been introduced. The first 
of them, Mo is zero everywhere except in the position fi$ of the diagonal. 
The second is the identity matrix except for the element /xo of the diagonal 
that is zero. Using now the trivial matrix identity 

:")P(!i)-P (13 » 

we can write (R Mo ) in terms of R 

(R Mo ) = ^ RS+ (14) 

by means of a "swapping" matrix that is inspired by the identity of 
Eq. [T3J The effect of acting to the left of the Bogoliubov amplitudes 
W, i.e. W/iq — WS/j, , is to swap the row fi of the U and V amplitudes. 
The structure of S^ is that of an identity matrix except for the rows and 
columns of the label no in both the U and V blocks. The simplifications 
implied by the introduction of 5 Mo can be seen for instance in the expression 
of the generalized density 

Tl^ = WR^W+ = W^ RW+. (15) 
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This tells us that the generalized density can be written in terms of the 
standard formulas (for instance p = VV T ), but using the new U and V 
matrices. More interesting is the variation of the energy at first order when 
the Bogoliubov amplitudes W are varied according to the most general 
canonical transformation (see Ref [24] for notation) 

W(Z) = W(0)e iZ . (16) 

where Z is an hermitian bipartite matrix 

/ yll 720 \ 

The variational parameters of the theory can be enumerated as: the com- 
plex matrix off-diagonal elements Z^ n with m > n; the diagonal Z^ m ; and 
the complex off-diagonal matrix elements Z^ n with m > n. The change in 
energy is given by 



SE = -Tr 2 [[R,H]Z] + 0(Z 2 ) (18) 



with 



and 



-W + {0)HW{0) = 



\-a* -(t + r) 

On the other hand, the change in energy for a blocked HFB state \4>fj, ) i s 
given by 

SE^ = ^Tr 2 [[R M0) H]Z] + 0(Z 2 ) (19) 

where we have replaced K by R Mo and with an H computed from the same 
density. Using the " swapping" matrix we obtain instead 

SE m = l -Tr 2 [[K,H M0 ]Z W ] + 0(Z 2 ) (20) 

with H, 1Q = S^MS+ = W+HWn and Z po = 5 M Z5+. For the Bogoli- 
ubov amplitudes, the following relation is helpful 

W(Z) lto =W(0) lto e a ^. (21) 
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In practical implementations of the gradient method the exponential in Eq. 
|2T1 is computed using the series expansion but corrected to have unitarity. 
A possibility is 

e iZ "o w (I + iZ^XI + Z^Z^) 1 / 2 . 

where the square root of the positive definite matrix is computed by means 
of the Cholesky decomposition. Others, based on Pade rational approxi- 
mations to the exponential have been explored |24) . 

The previous results are telling us that we can use exactly the same 
gradient formalism as in the even-even case but using as starting amplitudes 
W(0)fj, o . Obviously, the idea can be generalized to multiple quasiparticle 
excitations just by adding as many swapping matrices as quasiparticle 
excitations considered. 

These ideas are being extended to the expansion of the energy up to 
second order required for a "second order" (Newton like) gradient method 
and its descendants like the use of the inverse of the sum of quasiparticle 
energies + E v to damp the "high energy" components of the gradient 
G M „ as discussed in [34] . Although this little trick can not be used for fi- 
nite temperature systems ( K^ = is a necessary condition, not satisfied 
for finite temperature density matrices), work on an efficient implementa- 
tion of the gradient method using the inverse of two quasiparticle energies 
as a prc-conditioncr and valid for any situation (even-A, odd-A or finite 
temperature) systems is in progress ]33 . . 

4. Conclusions and perspective 

Although the standard BCS theory and its use in nuclear physics are both 
more than fifty years old, there are still technical issues, particularly re- 
lated to systems with an odd number of particles, that require further 
developments to simplify the systematic application of BCS/HFB and be- 
yond to nuclear systems all over the nuclide chart. In this contribution 
we have discussed two of them, one related to the overlaps of HFB wave 
functions required in theories beyond mean field and using the pfaffian of 
skew-symmetric matrices. The other focused on the gradient method with 
blocked HFB wave functions. In the near future, we hope to extend the 
pfaffian technique to finite temperature systems and make use of it to sim- 
plify the appliction of of the generalized Wick theorem. Also, approximate 
second order gradient methods will be extended to odd-A and finite tem- 
perature systems. 
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